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To contrast different generators for flow equations for Hamiltonians and to discuss the depen- 
dence of physical quantities on unitarily equivalent, but effectively different initial Hamiltonians, a 
numerically solvable model is considered which is structurally similar to impurity models. By this 
we discuss the question of optimization for the first time. A general truncation scheme is established 
that produces good results for the Hamiltonian flow as well as for the operator flow. Nevertheless, it 
is also pointed out that a systematic and feasible scheme for the operator flow on the operator level 
is missing. For this, an explicit analysis of the operator flow is given for the first time. We observe 
that truncation of the series of the observable flow after the linear or bilinear terms does not yield 
satisfactory results for the entire parameter regime as - especially close to resonances - even high 
orders of the exact series expansion carry considerable weight. 

PACS numbers: 03.65.-w, 05.10.Cc, 33.80.-b 

I. INTRODUCTION 
A. Flow equations 

Eight years ago, Glazek and Wilson! and independently Wegner 2 introduced a new non-perturbative method to 
diagonalize, renormalizc or simplify a given Hamiltonian. Whereas in high energy physics the method is known as 
"similarity transformations" , the term "flow equations" has been established in the solid-state community. The idea 
is conceptually simple: Instead of diagonalizing the Hamiltonian of the system by a single unitary transformation, 
one performs a continuous sequence of infinitesimal unitary transformations and thus induces a flow on the system 
parameters. The procedure is not constrained to specific symmetries nor to certain parameter regimes - but is accessible 
to any system described by a Hamiltonian. Thus, the method has been successfully applied to various models of solid- 
state and nuclear physics. Examples are dissipative quantum system a 3 i 4 i 5 ,the electron-phonon problem^ or the 
Hubbard model^. For a recent review on the flow equation method see Ref. |9j. 

The main advantage of the method is its flexibility. This is similar to the numerical diagonalization of a given matrix: 
There arc many different possibilities to reach the goal. One is free to choose the basis in which the diagonalization 
is performed and within a given basis one is free to choose the concrete series of unitary transformations that finally 
diagonalizes the matrix. Depending on the basis and on the concrete series of unitary transformations, convergence 
may be good or poor, numerical errors may be small or large. 

Similarly, many different flow equations can be formulated to diagonalize or simplify a given Hamiltonian. Even 
though all different flow equations are equivalent and will eventually lead to the same result, matters change as soon as 
approximations are involved. Typically one needs to cut the hierarchy of newly generated interaction terms and then 
neglect operators, which are assumed to be irrelevant. Yet, there is no satisfactory definition for irrelevant operators 
within the flow equation approach. Whether or not a contribution is irrelevant depends on the initial Hamiltonian 
and on the goal on wants to reach. 

Usually, approximations were justified when certain sum rules, mostly stemming from the invariance of commutation 
relations during the unitary flow, hold exactly or at least asymptotically 5 ^. In addition, exact relations between 
static and dynamic properties - as the generalized Shiba relation in the case of the spin-boson modeliS - can serve as 
justification for prior approximations 5 -. A general consistency check lies in the explicit investigation of the flow of the 
neglected operators. 

So far, a detailed discussion on optimization of flow equations is missing. With this work, we want to start to fill 
this gap by addressing the following questions: E.g., any initial Hamiltonian H implicitly depends on a number of 
parameters H = H(?p,9, ...) where the parameters are associated to certain unitary transformations Ui(ip), £/ 2 (#), ••• 
. Can the parameters ip, 9, ... be chosen such that a given flow equation scheme yields optimal results? A second 
variation that will be discussed lies in the arbitrary definition of the "diagonal" Hamiltonian, Hq - as mentioned above. 
Will different Hq yield similar results for physical quantities and is there an optimal Hq for all physical quantities - 
or does the optimal Hq depend on the physical quantity under scrutiny? 

Another fundamental question associated with the flow equation approach is connected to the observable flow and 
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has not been discussed in depth yet, either. For this, we note that in order to take advantage of the simple structure 
of the fixed point Hamiltonian, the observable has to be transformed as well - by the same sequence of unitary 
transformations that diagonalized the Hamiltonian. Since usually the continuous transformation is designed such 
that the diagonalization of the Hamiltonian is "optimal" , the observable flow is more likely to suffer from uncontrolled 
approximations. We will address the question if there is a scheme that optimizes both, Hamiltonian and observable 
flow and also compare the observable flow on the operator level. 

To do so, we will not proceed systematically but we will address these questions more specifically. Namely, we will 
consider an explicit model which is structurally similar to dissipative impurity models - but still exactly solvable via 
numerical diagonalization. We will call this model the Rabi model and it is presented in the next subsection. We use 
this model to test different approximation schemes and compare the results with the numerically exact solution. This 
strategy was first pursued by Richtcr in his diploma workii. We extend his work in various directions. One point is 
to investigate Hamiltonians where the reflection symmetry is broken. This is important if one wants to understand 
the mechanism of phase transition as being observed in the spin-boson modeU£. 

In Sec. II, we develop a general truncation scheme which yields good results over a wide range of the parameter 
space. Furthermore, we present a particular truncation scheme which leaves the Hamiltonian form-invariant during the 
flow. The question of the invariance of the flow equations with respect to the particular choice of initial Hamiltonians, 
provided that they only vary by a unitary transformation, is discussed. As a criterion for the quality of the flow 
equations, we look at the ground-state energy as function of the bias as an example for the flow of a parameter of the 
Hamiltonian. In Sec. Ill, we give a thorough discussion about the flow of observables. As reference we will not only 
investigate the expectation value of observables, but also compare the flow equation result with the exact solution 
on the operator level for the first time. For this, an expansion of the operator into a basis of normal ordered bosonic 
operators is given. In Sec. IV, we close with general remarks and conclusions. 



B. The Rabi Hamiltonian 



The specific model we use for our discussion of various realizations of flow equations and various approximation 
schemes is the spin-boson Hamiltonian with only one mode, which we will call the Rabi model in order to distinguish 
it from the spin-boson model with an arbitrary number of modes. The Hamiltonian is given by 

H = -^a x + ^<r s +uotfb + <r x ^{b + b1)+E . (1) 

Here, b^) denotes the bosonic degree of freedom and o~i with i = x,y,z are the Pauli spin matrices. They obey the 
canonical commutation relation [b, 6'] = 1 and the spin-1/2 algebra [cr^cr,-] = 2ieijkO~k- Since there is only one mode 
present, a numerical diagonalization is feasible by truncating the bosonic Hilbert space after n bosonic excitations 
with some fixed value of n. 

The model was first introduced in the context of spontaneous emission and absorption of atoms and due to its long 
history there exists an enormous amount of work that has already been published on this model. It is impossible to 
review or cite all these papers - a good overview may be found in the paper by Graham et aiii The model has also 
been discussed in connection with quantum chaosi^i^ and extensions of it can serve for the description of optical 
phonons interacting with two-level systems or quantum dots within a solid-state matris^. In the context of flow 
equations, the model has been discussed by Mielke- 1 -- using a set of flow equations that preserve the banded structure 
of the Hamiltonian. In the present work we focus on low energy properties of the Rabi model. Ref. is m some 
sense complementary to the present work since there the high energy modes were discussed. 

In the present work we are interested in general properties of the flow equation method. The reasons for us to 
investigate the Rabi model lie in the fact that it couples a two-level system to a "bath" resembled by the bosonic degree 
of freedom. And since we only consider one mode, the system is still exactly solvable via numerical diagonalization. 



II. FLOW OF THE HAMILTONIAN 



A. Setting up the basis 

In order to diagonalize the Hamiltonian |T|. we will perform a continuous unitary transformation. The flow 
equations are generated by the anti-hermitian operator rj which is canonically given by r\ = [Ho, H], where Hq defines 
the diagonal Hamiltonian 2 . Different choices of rj are possible as well. The flow equations are of the form 
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where both H and r\ depend on the flow parameter £. The choice 77 = [Ho , H] is likely to decouple the fermionic 
system from the bosonic system, and the fixed point Hamiltonian H(£ = 00) is then basically given by Hq where 
the asterisk indicates that the parameters of the initial diagonal Hamiltonian are in general renormalized. For a brief 
introduction, we refer to Appendix^ which treats the Rabi model with A = and motivates the approach given 
here. 

Obviously, different choices for Hq can lead to different flow equations. Another ambiguity stems from the fact that 
the initial Hamiltonian may differ by a unitary transformation. If we restrict ourselves to orthogonal transformations 
in the two-dimensional Hilbert space and to simple translations in the bosonic Hilbert space, i.e. 

the general initial Hamiltonian with respect to the Hamiltonian Q is given by 

H = -^a x + ^a z +u tfb+^(b + tf)+a x ^(b + tf)+a z ^-{b + tf)+E' , (4) 

where we introduced the following parameters: 

A' = Ao cos if> + e sin ip , e' = e cos ifi — Aq sin if) (5) 
A e = 8X0 , X x — — sinV'Ao , X z = cost/>Ao (6) 

E o = E + 2 ^ , ~e = e Q + 9^ (7) 

As was mentioned above, different generators and different unitarily equivalent Hamiltonians will lead to the same 
physical results if no approximations are involved. But the above model is not solvable analytically and therefore 
approximations become necessary. In the special case of the Rabi model the flow equations will generate an infinite 
series of new coupling terms which cannot be summed up formally to yield a closed expression^. 

In this section we will first only take coupling terms into account which are linear in the bosonic operators and 
have real coefficients. This means that with respect to the initial Hamiltonian Q), only the term ia y (b — w) will be 
newly generated which resembles the lowest order of the polaron transformation (see e.g. Ref. Il2t) . Using a generator 
which is not of the simple form 77 = [Ho,H], we will also discuss flow equations which leave the initial Hamiltonian 
form- invariant. We are able to show analytically that the fixed point of the flow equation is independent with respect 
to the (distinguished) unitary transformation. 



B. Flow equations with respect to the Canonical Generator 



In the following subsection we will discuss flow equations which are obtained by employing the canonical generator 
7/ = [Hq, H]. This gives rise to new interaction terms. The truncated Hamiltonian shall be given by 

A f \ e \ x \y \ z 

H=--a x +-a z +u tfb+—(b + tf)+a x — (b+tf)+ia y — (b-tf)+(j z — (b + tf)+E , (8) 

where all parameters but the bath energy loq are explicitly ^-dependent. The above Hamiltonian represents the most 
general Hermitian operator which includes all possible interaction terms acting on the underlying Hilbert space up to 
linear bosonic operators with real coefficients. 
The flow shall be governed by the generator 

77 = layri ^ + 77 e (6 - b r ) + a x T] x (b - tf) + ia yr] y (b + tf) + a zV z (b - b r ) (9) 
= fi°' v + rf + 7f + 7f + rf , 

where the parameters 77 ' 1 ', rj e , r\ x , 77^, and jf are ^-dependent and will be specified later. 

The above generator represents the most general anti-Hcrmitian operator which includes all possible operators 
acting on the underlying Hilbert space up to linear bosonic operators with real coefficients. 
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1. Setting up the flow equations 

The commutator [77, H] yields the following contributions: 

[fj°' v , H] = -a z Arf' v - a x erf' v + a z rf> v \ x (b + &t) - a x rf' v X z {b + (10) 
[rf , H] = 7] e u; Q (b + fct) + 7 fy + ax rf\* + a z jfX z (11) 
[if, H] = -ia y er) x (b - 6+) + a x r, x w {b + tf) 

+ o x rfX e + rfX x - a z rfX y {b - &t) a - ia yV x ^-{(b - b+), (6 + fe f )} (12) 

[rf , H] = -a z Ar) y (b + 6*) - a x a?{b + 6*) + ia y r/ v Lu (b - 6 f ) 

+ a z jf A x (6 + & f ) 2 + ?f A y - <7 X ^A Z (6 + 6 f ) 2 (13) 
[jf , H] = -i^A^i - 6 f ) + atfwoQ) + & f ) + a^X" 

+ *'<V7* y {(& - & f ), (6 + & f )} + ^ Z A^(& - &t) a + ?f A z (14) 

{.,.} denotes the anti-commutator. As can be seen in (112(1 - l|14jl. the flow equations generate terms which are 
bilinear in the bosonic operators and we will need to find a suitable procedure how to include these terms in the flow. 
Kehrein, Mielke, and Neu^ proposed to neglect these terms after normal ordering them with respect to a bilinear 
bosonic Hamiltonian. Since we allow the initial Hamiltonian to differ by a shift in the bosonic operators, we need to 
include this generalization also in the normal ordering procedure, i.e. we will normal order with respect to the shifted 
bosonic mode 

£=&+- , (15) 

with the linear shift 6 to be determined later. To close the flow equations we will thus neglect the normal ordered 
operators 

O x = -a x rf>X z :{l + Pf : , 2 = a z rfX x : {b + Pf : , (16) 
O a = <r x ri z Xy ; (5- V) 2 : , O4 = -o z rfX y : (b - V>) 2 : , (17) 

o, = l a y (v z Y ~ ^ X y) : " & f )> (* + 5t )> : • ( 18 ) 

Normal ordering is now defined as : (b + P) 2 := (b + & f ) 2 - 1„, with 1„ = ((& + P) 2 ) = 1 + 2n, and n = (e 0UJ ° - 1) _1 
being the Bose factor. Notice that the temperature enters in the Hamiltonian flow through normal ordering. In the 
following we will only consider T = 0, i.e. 1„ = 1, but we will nevertheless keep track of this distinction. 

Like in the case of flow equations for impurity systems^, the above truncation scheme has the effect that the bosonic 
energy loq is not being renormalized during the flow. 

With =§■ = [77, H] we obtain the following flow equations: 

d t A = 2er)°' y - 2r/ e X x - 2rfX e + 2(r) z X y + ri y X z )l n - 2r/ y X z S 2 (19) 

d e e = -2Ar 1 ' y + 2r, z X e + 2(r l y X x +r 1 x X y )l n + 2r, e X z -2ri y X x 6 2 , d e X e = 2 V e uj (20) 

de\ x = —2erj y + 2r/ x Lu — 2rf' y X z + Arj y X z S , d t X y = -2Aif - 2erf + 2r) y u - 2rj z X x S + 2r/ x X z 6 (21) 

diX z = -2Arj v + 2ti z ujo + 2rf> y X x - 4ri y X x S , dtE = rj e X e + rj x X x + rfX y + T) Z X Z (22) 

With A e = rf = 0, an obvious invariant is given by Inv = A 2 + e 2 + X x2 + X y2 + X z2 — AEujq. To investigate the flow 
equations further, one has to specify the constants and initial conditions. To do so we will choose different diagonal 
Hamiltonians Hq, and we will contrast the resulting flow equations by means of the ground-state energy of the system. 



2. Determining the Canonical Generator 

An obvious choice for the diagonal Hamiltonian is given by Hq = — —a x +Uob^b. The canonical generator rj = [Hq, H] 
is of the form © with ?; ^ = Ae/2, if = -u X e /2, rj x = ~uj X x /2, rf 1 = (AA z -w A y )/2 and rf = (-oj X z + AX v )/2. 
We will refer to the flow equations with this particular choice of rj as Version a. 
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Another choice for the diagonal Hamiltonian is given by Hq = j<j z +u)otfb. The canonical generator 77 = [Hq, H] is 
of the form © with r)°' y = -Ae/2, r/ e = -w A e /2, 77* = (-u \ x + e\ v )/2, rf = (e\ x - uj W)/2 and rf = -lu X z /2. 
We will refer to the flow equations with this particular choice of i] as Version b. 

The third choice for the generator which we will investigate in the following combines the two previous choices, i.e 
r\ = [Hq, H] with Hq = —^u x + ^az+LOob^b. The canonical generator r\ is of the form © with r)°' v — 0, rf = — uj$\ e /2, 
rf = (-woA* + eA J/ )/2, rf = (A\ z + eA* - ui \ y )/2 and rf = (-u \ z + AA 3/ )/2. We will refer to the flow equations 
with this particular choice of r\ as Version c. 

There are other possibilities for the diagonal Hamiltonian which include coupling terms. We could e.g. choose 
H a =u tfb + a z f(b + 6+), since this Hamiltonian is also exactly solvable, see Appendix^] Another possibility is to 
choose the Jaynes-Cummings Hamiltonian as Hgl& (see also Appendix [D| , which was done by Richter—. In this work 
though, we want to confine ourselves to the versions given above. 



3. Determining the Bosonic Shift 



We now want to determine the newly introduced bosonic shift 8. The procedure is not unambiguous, but we are 
led by formally diagonalizing the Hamiltonian as follows: 

A' e' „+ v X3 xv w, ^ A* , w 

= ~~2 Gx + 2 a * + ^ +2^ a ^ + lay 2^ ){h + 2^ a *2^- 1<T *2^ + E ' (23) 

3 j 

with A' = A + - E>L e ' = e- - and E 1 = E - K^- - and summation is over 7 =e,x,z with 

cr e = 1. Decoupling the fcrmionic and bosonic Hilbert space, we thus obtain the ^-dependent shift 

* = £fo>- ■ ( 24 ) 

The fermionic expectation values can be evaluated directly with respect to the effective Hamiltonian H v = —^-<J X + 
\o z to yield 

(a x ) = A'/R' , (a z ) = -e'/R' , with R' 2 = A' 2 + e' 2 . (25) 



There is also a self-consistent possibility to determine the system expectation values. For that we will formulated the 
Hamiltonian with respect to the shifted mode b = b + S/2. The renormalized "one-particle" parameters are then given 
by 

A=A + \ X 5 , e = e-X z 5 . (26) 

Evaluating the system parameters now with respect to the system Hamiltonian H p = ——o x + \oz and still assuming 
the bosonic shift as given in Eq. (|24p. we obtain the following self-consistent equations: 

(a x )=A/R , (a z )=-e/R , with R 2 = A 2 + 1 2 . (27) 

In this work, we will restrict our investigation to the bosonic shift of (|24|) and to these two procedures of determining 
the fermionic expectation values. But there are other possibilities of evaluating the bosonic shift or the expectation 
values. One way is e.g. to couple the flow of the system parameters with the flow of the observable by imposing 
that a certain sum rule holds exactly (see next section). This condition will determine the bosonic shift. In the next 
section we show that the sum rule for the x- and z-component of the Pauli matrices is quadratic in the bosonic shift. 
But since we restricted ourselves to real shifts, there might be no solution. Even if we allowed imaginary coefficients 
in the evolution of the Hamiltonian, a solution would not be guaranteed since the sum rule would then relate the 
complex shift 6 with its complex conjugate 5*. Numerical investigations indicated that the bosonic shift 6 cannot be 
chosen such that a certain sum rule holds exactly. It is left open, how this effects the stability and reliability of the 
flow equation approach. 

Finally, we want to point out that the procedure of determining the expectation values can significantly alter the 
behavior of the flow equations. In case of the spin-boson model it is shown^S, that an infinitesimal bias resembles a 
relevant perturbation, i.e. die cx e for small £, if one chooses the expectation values directly whereas it resembles a 
irrelevant perturbation (die cx — e) if one chooses the self-consistent scheme. 
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FIG. 1: The ground-state energy E^ E obtained by different canonical generators with tp = and 9 = for Ao/loq = 0.5 (left 
hand side) and Ao/^o = 1.5 with Ao/u>o = 1 as a function of the bias eo relative to the exact ground-state energy Eg X , shown 
in the panel. The bosonic shift is set zero throughout the flow, i.e. 5 = for all £. 



4- Numerical Results 



We want to analyze the quality of the above flow equations by means of the ground-state energy E g of the system 
as a function of the external bias eo- These results are compared with the numerically exact solution obtained via 
numerical diagonalization. Since the bosonic mode is left un-renormalized, the energy scale is given by loq. For the 
coupling constant we choose Ao = loq, i.e. we are not in the perturbative regime. 

We will first consider the flow of the initial Hamiltonian with 9 = and ■0 = 0. We will also set S = for all t. 
In Fig. ^the ground-state energies Eg E obtained from the different canonical generators are shown. Calculations 
are done for two different tunnel-matrix elements Ao/^o = 0.5 (left hand side) and Aq/ujo = 1.5 (right hand side), 
the first below and the second above resonance. Resonance in the un-perturbed system is defined by Ao/wo = 1- All 
results are in good agreement with the numerically exact solution. Still, differences occur in the non-trivial regime 
where the bias eo is below or around the energy scale given by w . In the panels, the exact ground-state energies Eg X 
are displayed. 

We now turn to the flow equations obtained by employing the generalized normal ordering procedure, i.e. we set 
6 = J^j (crj) X' /u>o- The results for the different generators are shown in Fig. [3 The expectation values are determined 
directly according to Eqs. 1)2 5|) (left hand side) and self-consistently according to Eqs. (|27|) (right hand side). There 
is a systematic improvement to the results of Fig. ^ where S was set zero for all £. The best results are obtained by 
the generator of Version b and determining the expectation values self-consistently. 

Finally, we want to investigate the dependence of the flow equations on the unitarily equivalent, but different 
representations of the initial Hamiltonian, labeled by ip and 9. For this we choose the generator of Version b and the 
bosonic shift of l|24|) with the direct evaluation of the expectation values according to Eqs. (|25|) . On the left hand 
side of Fig. [3]we vary ip with 9 = 0; on the right hand side of Fig. [3]we vary 9 with if) — 0. 

As can be seen, there are differences with respect to the initial Hamiltonian. For i/j = 7r/4, there is a big deviation 
from the exact value in a small region around eo ~ 1.5 with a maximum of 1.2. In this region the fixed point 
Hamiltonian H(£ = oo) varies from the "normal" fixed point Hamiltonian and the ground-state energy is mostly 
determined by E(l = oo). This is also the case for 9 < — 1 (not shown) where the regions of large deviations depend 
on 9. Still, we observe a certain invariance with respect to the initial Hamiltonian keeping the crude truncation scheme 
in mind. 
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FIG. 2: The ground-state energy Eg E obtained by different canonical generators with ip = and 6 = for Ao/cjo = 1.5 with 
Ao /wo = 1 as function of the bias eo relative to the exact ground-state energy Eg X , shown in the panel. The expectation values 
for the bosonic shift 8 = ~Y2j{<?j)\ 3 /loq were evaluated directly according to Eqs. 12511 (left hand side) and self consistently 
according to Eqs. 1271 (right hand side). 



From the considered parameters, the best results are obtained for ip = and 9 = —0.5. Of course, it would be 
desirable to give an objective scheme how to choose the representation of the initial Hamiltonian that yields the best 
result for the ground-state energy. This had to be left open. 



C. Flow equations with respect to a Form-Invariant Flow 

As was mentioned above, the canonical generator r\ = [Hq,H], in general, gives rise to new interaction terms. In 
order to avoid this complication, Kehrein, Mielke, and Neu pursued a different strategy to set up the flow equations, 
namely they chose the generator rj such that the Hamiltonian remains form- invariant. To assure that the initial 
Hamiltonian of Eq. (QJ remains form-invariant, we set 5 = 0, and the constants of the generator of Eq. [5] have to 
satisfy the following relations: 

jf = , -etf + ri x uj Q - rf> v \ z = , -Arf - erf + rfuj a = (28) 

This guarantees that A e , A^ and X v are not being generated. With these relations, the parameters are defined up to a 
common factor /. If one chooses jf = —oJo\ z f/2, one finds rj°' v = eAf/2, if = 0, rf = and r\ v = ~AX Z f /2. With 
this choice, all neglected operators except of 0\ vanish. One obtains the following coupled differential equations: 

a,A = -AA- 2 /ln + Ae 2 / , d e e = -eA 2 f 
d e \ z = X z (A 2 -u J 2 )f , d e E = -io \ z2 .f/2 

For the numerical calculations, we set / = 1 and refer to this set of flow equations as Version d. 

We want to consider the form-invariant flow after having performed a unitary transformation on the two-dimensional 
Hilbert space which diagonalizes H p = —^Tx + ^&z — ► -f cr z with R 2 = Aq + e\. This is achieved by choosing 
tani/; = — Ao/eo- If we thus want to avoid the generation of A, A e , and X y as defined in J8J, we set 5 = 0, and the 
parameters of the generator have to satisfy the following conditions: 

Rrf> v + jf\ z \ n = , 7f = , -Ar] z -Rr] x +7ftj = (30) 
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FIG. 3: The ground-state energy Ef E obtained by the canonical generator of Version b for Aq/cjo = 1-5 and Xq/ljo = 1 as 
a function of the bias eo relative to the exact ground-state energy Eg X , shown in the panel. The expectation values for the 
bosonic shift 5 = X^j( <J j)-^ J / LJ o were evaluated directly according to Eqs. 12511 . The initial values were 6 = and various ip 
(left hand side) and tp = and various 8 (right hand side). 



Again the parameters of the generator are only defined up to a common factor. Choosing rf — 0, rf = —LUoX x f/2, 
rf = —uj X z f/2 renders O5 zero and yields rf' y = X x X z fl n /2 and rf — —RX x f/2. Thus all neglected operators but 
0\ and 02 are zero. We obtain the following How equations: 

d t R = -RX x2 fl n , d(E = -u! (X x2 + X z2 )f/2 (31) 

d e X x = -u 2 X x f + R 2 X x f - X* 2 X x fl n , d e X z = -u 2 X z f + X x2 X z fl n 

The set of equations in (|31|l is equivalent to the set of equations in (|29|l . This can be seen by introducing "new" 
variables A' = X X R/X' , e' — X Z R/X' and A' = X x2 + X z2 and setting up their differential equations, which coincide 
with l|29|) . This demonstrates that keeping the Hamiltonian form-invariant during the flow preserves the unitary 
equivalence with respect to the initial Hamiltonian for this special unitary transformation. 

If we want the initial Hamiltonian to remain form-invariant during the flow after having shifted the bosonic mode 
by 9, the constants have to satisfy the following relations: 

-erf +t] x lj -t] ' v X z + 2rjyX z 5 = (32) 
-Arf - er) x + ^luq + 2ij x X z 5 = (33) 

After the shift, X e is naturally generated which was not present in the previous schemes. In order to compare the 
flow equations with the above versions, we have to couple the flow of A e with the flow of X z , i.e. A e = 8X Z . This sets 
another condition on the parameters of the generator, i.e. if = —9Ari v /ujo + 9r] z . If we further choose rj y = —AX z f/2 
we obtain rf = -u \ z f/2, rf = 0, rf' y = eAf/2 - AX z 6f and if = 8A 2 X z f/(2uj ) - 9uj X z f/2 with the factor / to 
be determined later. With e = e — 9X z2 /wo, this yields the following flow equations : 

d i A = ~AX z2 fl n + A(e + X z (5-9X z /u J o)) 2 f , d e e = ~A 2 ef + 2A 2 X Z (5 - 9X z /u )f , 
d e X z = X z (A 2 -u, 2 )f , d e E = ~cj X z2 f/2 + 9 2 (A 2 -u, 2 )X z2 f/(2tj a ) 

Recalling the initial condition of the energy shift E' Q = Eo + 9 2 X^ / (Aujq) defined in Eq. Q we see that the flow 
equations are equivalent to the flow equations of Version d if we set / = 1 and 8 = X e /uiQ = 9X Z /ujq. This choice of 
the ^-dependent shift coincides with the expression l|24|) if we set (a z ) = 0. 
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FIG. 4: The ground-state energy E^ E obtained from the form-invariant flow for Ao/u>o = 0.5 (left hand side) and Ao/ljo = 1.5 
with Ao/wo = 1 as a function of the bias eo relative to the exact ground-state energy Eg X , shown in the panel. The primed 
versions are including the flow of the neglected operators 0\ and O2 (see text). 



This is a remarkable result. It shows that if one imposes invariance of the flow equations with respect to unitarily 
equivalent initial Hamiltonians and chooses the truncation scheme that leaves the Hamiltonian form- invariant, the 
flow equations are uniquely determined. It also shows that normal ordering with respect to the ^-dependent mode 
b = b + S leads to reasonable results. 

We now want to check the quality of the form-invariant truncation scheme. In Figure 01 the ground-state energy 
Eg E obtained by the set of equations lf^9l is shown relative to the exact ground-state energy Eg X for two different 
tunnel-matrix elements Aq/loq = 0.5 (left hand side) and Ao /u)q = 1.5 (right hand side). Drastic deviations from 
the exact result are seen in the regime eo/^o > 1. This means that the neglected operator 0\ of Eq. (|l(jfl becomes 
relevant and has to be taken into account. 

In order to demonstrate that the flow equations can be improved systematically, we will now consider higher 
order terms of the bosonic operators in their normal ordered representation. For the normal ordering procedure see 
Appendix IE! Since we set S = for all £, normal ordering is defined with the respect to the unshifted mode, i.e. b = b. 
Redefining 0\ = <j x k\ : (b + b^) 2 :, the commutator [77, Oi] yields 

[77, d] = 2a z rf fi Ki : (6 + 6 f ) 2 : +2a z r) y K 1 (: (b + 6 + ) 3 : +2{{b + b^f) : (b + tf) :) 

+ 2ia y r) z Ki : {b - 6 f )(6 + tff : . (35) 

We first neglect the trilinear operators and the bilinear operator of type O2 (see Eq. (|16fl ). The extended flow 
equations then read (/ = 1) 

d e A = -AA z2 l„ + Ae 2 , d e e = -eA 2 , dim = AA z2 /2 

d e X z = \ Z (A 2 - u)l) +4A z A Kl l„ , d e E = -uj \ z2 f/2 . " A " 

We will refer to this set of flow equations as Version d'. 

To see if this improvement is systematic we will now include also the corrections that come from the neglected 
operator of type 02- Redefining C 2 = cr z K2 : (b + b^) 2 :, we obtain similar commutator relations for [77,02] as we got 
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in 123): 

[q, 2 ] = -2a x rf •% : (b + tf) 2 : -2a x r, y ^{- (b + & f ) 3 : +2((6 + fc*) 2 ) : (b + tf) :) 

-2ia y r) z K 2 : (fo- 6+)(fo + fot)2 . ( 37 ) 

The effect of including the operator 2 in the flow equations is the following: The conditions for the constants of the 
generator that assure the form-invariance of the Hamiltonian slightly change, see Eq. Q28JI. The flow equations thus 
read (/ = 1) 

d t A = -AA 22 1„ + Ae(e + 4 K2 ) , die = — (e + 4k 2 )A 2 

d t \ z = A Z (A 2 - uol) +4A 2 Aki1„ , d t E = ~uj q X z2 f /2 (38) 
= AA z2 /2- A(e + 4K 2 )K 2 , cW = A(e + 4k 2 )ki . 

We will refer to this set of flow equations as Version d". 

In Fig. 0]one sees that the extended flow equations yield a systematic improvement ranging over the whole parameter 
space. Nevertheless, the agreement with the exact result remains rather poor for eo/wo > 1. Only if one considers the 
renormalization of the bath mode u>q, one obtains results within a few percent relative error over the whole parameter 
range. Regarding the spin-boson model, it is preferable to employ the canonical generator since the bath modes 
remain unrenormalized in the thermodynamic limit^i. 

III. FLOW OF OBSERVABLES 

We will now investigate the flow of observables. In order to characterize the quality of the flow equations, normally 
sum rules are derived expressing the fact that of = 1 or (anti-)commutation relations should hold for all £ with 
i = x,y, ,a£i As will be pointed out in the end of this section, these sum rules can be misleading. We will therefore 
contrast the expectation value (a z ) as it follows from the flow equation approach with the numerically exact solution. 
Furthermore, we will compare the flow equation results with the numerically exact fixed point of the operator flow on 
the operator level. To do so, we will give a unique decomposition of the fixed point operator into a basis of normal 
ordered bosonic operators. 

A. Flow Equations for the Pauli Matrices 

In order to take advantage of the simple form of the fixed point Hamiltonian when calculating expectation values of 
observables, the observable has to be subjected to the same sequence of unitary transformations as the Hamiltonian. 
The flow equations for the Pauli spin matrices thus read didi = [77, ai\. Again the flow equations generate an infinite 
series of operators and one needs a suitable truncation and decoupling scheme. The i-component of the Pauli spin 
matrices as a function of the flow parameter £ shall be given by 

ai{t) = 9i {i)a x + h(£)a z + fi(£) + a xX x ' l (£)(b + 6 f ) + in vX v ' 1 (£)(b - 6+) + a zX ^(£)(b + 6 f ) , (39) 

with i = x, z. We want to emphasizes that the constant term /$ is indeed generated even though it seems to contradict 
the theorem of the invariance of the trace under unitary transformations. A short discussion is given in Appendix iBl 
The flow of the y-component of the Pauli spin matrices is given by 

i*v(jt) =g v (i)ia y + o x X x ' v {b-tf)+<j zX z >y(b-tf) . (40) 

These are the most general expansions up to linear bosonic operators with real coefficients that can evolve from the 
Pauli spin matrices under the flow equations, i.e. from oi(£ = 0) = (7$. 
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The commutator [77, Gi\ with i = x,z yields the following contributions: 

[rj°-y, ^{1)] = 2a z9l T 1 ° y - 2a x h iV ' y + 2rr ; //' " x " ( + 6+) 

-2<7 x f/ ) 'V i (& + & + ) 
{r ] e ,a l (£)}=2<T x r 1 e X ^ + 2a z r 1 e X ^ 

[V x ,<n(l)} = -2ia y h iV x (b - 6t) + 2^ X ^ 4 - 2a x rfx v '\b - & f ) 2 
- i(7 ^ x ^{(6-6t),(6 + 6t )} 

[r) y , <Ti(l)] = 2a z9i r) y {b + &*) - 2h i( T x r) y {b + 6+) 

+ 2a z rf X x '\b + & f ) 2 + 2r?V < - 2tr x » 7 V < (& + & f ) 2 

[Tf,vM = 2iVy9iV z (b - & f ) + ^^{(b - & f ), (6 + & f )} 
+ 2 ( r ;c ^ X ^(&-& t ) 2 + 2^x z ' 1 

The commutator [77, icr^] is given by: 

[J? * <7„(*)] = 2a^°^ x ^(6 - 6+) - 2a,,, :l 'V "(>, - &t) 
[77*, a y {£)} = -2u z9y rf{b - 6+) - 2icT v T, x X z ' y (b - 6*) 2 

<r„ (*)] = a.rf X ^{(& + fc f ), (6 - b r )} - a xV y X z ' y {(b + (6 - & f )} 
[77* , a y (£)} = 2a x g v f] z (b - &t) + 2i<7 J/ ^ X x ' 3/ (6 - & f ) 2 



(41) 

(42) 
(43) 

(44) 

(45) 



(46) 
(47) 
(48) 
(49) 



Again, {., .} denotes the anti-commutator. To understand which operators can transform into one another, we give 
a list of operators and their behavior under parity transformation (P) and Hermitian conjugation (H) (x = (b + tf), 
p=(6-fct)) : 





1 


Ox 


ICfy 




X 


p 


a x x 




lUyX 


iCFyP 


a z x 




p 

H 


+ 
+ 


+ 
+ 




+ 


+ 
+ 




+ 




+ 


+ 

+ 


+ 

+ 


+ 



In order to close the flow equations, we neglect normal ordered bosonic bilincars where normal ordering is defined 
with respect to the shifted bosonic mode b = b+S/2. Thus, one obtains the following set of linear differential equations 



for the z-component of the Pauli spin matrices with i = x,z: 

d l9l = -2hiT] ' y + 2rf X x - 1 - 2r) z X y> l \ n - 2rf X z ^\ n + 2f 1 y X z < I 6 2 (50) 

d t hi = 2 9l if y + 2rf X x n n + 2 V x X y ' i l n + 2,7V 4 - 2if x ^8 2 (51) 

defi = 2r ) x X x ' i + 2r 1 y x v ' 1 + 2r 1 z x z > t (52) 

d eX x ' 1 = -2hiT]y - 2rf > ' y X z ' i + Aafx*^ (53) 

d eX y,i = 2 gi if - 2hitf + 2rf X z H - Tsf-fH (54) 

d lX z ^ = 2 9lV y + 2^y x ^ - ±if x x H (55) 

The flow equations for the y-component read: 

d l9y = -2 V z X x -yi n + 2r, x X ^l n (56) 

d lX x -y = 2 9y r) z - 2r] ' y X z ' y - 2rj y x z ' y S , d eX z ' y = -2 9y r,^ y x + 2r]°' y x x ' y - 2ri y X x ' y S (57) 



If no approximation was made, = 1 would hold for all £ and i = x,y,z. Taking the expectation value with 

respect to the bilinear Hamiltonian of the shifted modes the relation should hold approximately for i = x, z: 

= gf + hj + ff + w + x y ^ + x z ' l x z - l )i n 

+ 2{ 9i {a x ) + hi{a z ))fi + 2{ X x ' i {a z ) - x*>*»X^ (58) 
+ (X X - 1 X X ' 1 + X Z - Z X Z '1S 2 - 2((g+(a x )f) X x '* + (h + (cr z )f) X z ^8 
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FIG. 5: The expectation value {a z ) FE obtained by different canonical generators with i/> = and = for Aq/cjo = 1-5 with 
Ao/wo = 1 as function of the bias to relative to the exact expectation value {a z ) ex , shown in the panel. The expectation values 
for the bosonic shift S = X^j( <7 j)^ : '/ LJ o were evaluated directly according to Eqs. 125H (left hand side) and self-consistently 
according to Eqs. 1271 (left hand side). 



For the y-component we obtain: 

(^W> = 9l + (X X ' V X X ' V + X z ' v X z ' y )ln « 1 (59) 

Other conservation relations follow e.g. from the commutator \a x {t),a z (t)] = —2ia y (£). These relations can be used 
to assess the validity and the quality of the flow equations but they cannot assure whether the scheme will yield the 
correct results. We will comment on this point at the end of this section. 



B. Numerical Results for the Expectation Value of a z 



Measurable quantities other than the ground-state energy are determined by means of the operator flow. In this 
subsection we will discuss the expectation value (a z ) as it follows from the different versions of the flow equation 
approach. The expression is given by 

(o*> = *Wz{l = oo)>* = g(£ = oo)*(o- x )* + h(t = oo)*(<r z )* + f(£ = oo) . (60) 

Here, *(...)* denotes the ground-state expectation value with respect to the fixed point Hamiltonian H(£ — oo). 

In Fig. 0we contrast the results for the different generators which were discussed in the last section, (a z ) FE , with 
the numerically exact solution (a z ) ex . We choose tp = and 8 = for the initial Hamiltonian and we will employ the 
flow equations obtained by the generalized normal ordering procedure, i.e. S = Y^j^j)^ 3 V w 0' 

On the left hand side of Fig. the expectation values in the expression of S are determined directly according 
to Eqs. On the right hand side of Fig. the expectation values are evaluated self-consistently according to 
Eqs. [Z7\ For cq/ojo > 1, the best results are obtained by the generator of Version b with the direct evaluation of the 
expectation values. But deviations from the exact solution in the region eo/wo — 1 are significant. In the latter region 
the generator of Version c yields the best results. We recall that the ground-state energy was best approximated by 
the generator of Version b with the self-consistent evaluation of the expectation values entering the bosonic shift 8. 
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FIG. 6: The expectation value {a z ) FE obtained by the canonical generator of Version c with 8 — ~}2 (aj)X^ /uiq for Ao/wo = 1.5 
and Ao/i^o = 1 as a function of the bias eo relative to the exact expectation value {a z ) ex , shown in the panel. The parameters 
of the initial Hamiltonian are given by 9 = and various tp (left hand side) and tp = and various 9 (right hand side). 



This demonstrates that the "best" generator and "best" procedures of taking account of the neglected terms might 
depend on the physical quantity under consideration. 

We will now also include the initial unitary transformation on the two-dimensional spin-Hilbert space, label by tp 
and the initial bosonic shift 9 in our discussion. We will use the generator of Version c with the direct evaluation of 
the expectation values. On the left hand side of Fig. we vary tp with 9 — 0; on the right hand side of Fig. Ejwe 
vary 9 with ip = 0. 

Regardless the initial Hamiltonian, the flow equation results differ from the exact solution in the region eo/too < 2. 
But some initial Hamiltonians provoke more significant deviations than others. Good results over the whole parameter 
space are obtained by combining non-zero values of tp and 9 which "compensate" their errors, e.g. tp — ir/32 and 
= —0.2. Nevertheless, we were not able to given an objective procedure how to choose the optimal initial Hamiltonian 
- a priori. 



C. Operator Fixed point 

It is possible to compare the exact results with the flow equation approach not only on the spectral but also on the 
operator level. For this we have to diagonalize the Hamiltonian in this basis in which the corresponding "diagonal" 
Hamiltonian Ho of the flow equation approach is diagonal. Let now Hd = U HU* denote the diagonalized Hamiltonian, 
then a* = U<JiU' f is the operator to be compared with o~i(l — oo) stemming from the flow equation approach, with 
i = x,y, z. To do so we will decompose a* in a set of operators which are created by the corresponding flow equations. 

If one uses an expansion which is normal ordered in the bosonic operators the decomposition can be obtained 
numerically without any approximation^. The reason for this is that the bosonic ladder operators cannot compensate 
each other and then act on lower bosonic subspaces. To make this more explicit the general matrix structure of a 
normal ordered operator consisting of N bosonic operators is shown on the left hand side of Figure [7| taking the 
set {\u}} as basis with \v) = (tf) u /ViA\0) and 6|0) = 0, v being a positive integer. The dark area contains non-zero 
entries whereas the white area contains no entries. In case of a non-normal ordered operator the white, upper left 
triangle would also contain non-zero entries. 

As an explicit choice of the operator basis for real symmetric operators like o x and a z we choose the set {o : 
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FIG. 7: Left hand side: Matrix structure of a normal ordered bosonic operator consisting of N bosonic operators with respect 
to the canonical basis (see text). The dark area indicates non-zero entries. Right hand side: The dark area indicates the matrix 
elements of an arbitrary matrix which are uniquely determined by normal ordered operators consisting of up to a certain 
number of bosonic operators (see text). 



(b + tf) n (b - tf) 2m :,o' : (b + tf) n '{b - &t)2m'+i where Q = \^ x ^ z an d ' = ia y . The operator basis for real 
antisymmetric operators is obtained by interchanging o and d . In the following we will only consider the flow of real 
symmetric operators. The results also hold for the real antisymmetric case. 

We want to decompose a real symmetric operator into a set of finite operators. Considering all oper- 
ators of the basis given above with less or equal than 2iV-bosonic operators, we obtain a finite basis of 

3 J2m=o X)«=cT + J2m'=o Xm'=i~ m = + l)(4iV + 3) operators. Summing up the independent matrix elements 
which are uniquely determined by the normal ordered operators containing up to 2N bosonic modes, we obtain 
J2n=o 2 ( 4n +!) + ! = 4(iV + 1)N + 3(N + 1) = (N + l)(4N + 3). These independent matrix elements are located at 
the upper left triangle of the matrix, indicated as dark area on the right hand side of Figure [7| 

In order to complete the discussion we also consider all operators with less or equal than (2iV+ l)-bosonic operators. 

We then obtain a basis with 3 J2m=o Enio"™'" 1 " 1 +Em'=oEn'=(j' m = + l)(4iV + 7) operators. Summing up 
the independent matrix elements which are uniquely determined by the normal ordered operators containing up to 

2N + 1 bosonic modes, we obtain J2n=o 2 ( 4n + 3 ) + 1 = A ( N + V N + 7 ( N + I) = {N + 1)(4JV + 7). 

We thus obtain the same number of independent matrix elements and basis "vectors" . This confirms that our 
basis is complete and linearly independent as we take N — > oo. Secondly, this shows that the first (N + 1)(4JV + 3) 
coordinates of a real symmetric operator with respect to a finite basis of operators up to 2N bosonic operators are 
left unchanged if one goes over to a finite basis including 2N + M bosonic operators (M > 0). 

We can thus exactly determine the coefficients of our basis up to any number of bosonic excitations N which a* 
is composed of. This shows that choosing a set of normal ordered bosonic operators as a basis yields a systematic 
approximation of any operator. If one is only interested in the system dynamics at low energies it thus suffices to 
consider only up to N bosonic operators with N — 2 say. 

To determine the coefficients numerically one has to work with a specific basis. Up to now we have only specified 
the basis of the bosonic Hilbert space. Choosing Hq = —^-a x + LVob^b to be diagonal we are led to the basis {|e, v)} 
with the first quantum number e = 0, 1 denoting the eigenstates of a x and the second quantum number denoting the 
eigenstates of tfb. Choosing Ho — |cr 2 + Lo^b (Version 16) or the diagonalized representation Hq = ^a z + ui^b^b of 
Version lc, we would choose the first quantum number e = 0, 1 to denote the eigenstates of a z . 

Considering all operators with less or equal than 2A r -bosonic operators, we end up to solve a linear equation Ax = b, 
with A being a quadratic matrix and x,b being vectors with dimensions (N + 1)(4JV + 3). The coefficients of the 
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matrix A are obtained by the following matrix representations of normal ordered bosonic operators: 

n / v. 2m ,„ v 

(e lM |o : (6 + 6t)n (6 _ 6 t )9 m : | e » = ( e | |e')£ (?) £ , ) (-1) 2 



fc=0 v 7 i=o x 



\2m-l 



x 8(1/ - k - I) J {N k l){ J {N k , <S M ,„ +n+2(m - fe -z) 



(61) 



The vector & on the right hand side of the linear equation is given by the (N + l)(4A^+3) independent matrix elements, 
located at the dark area of the matrix of the right hand side of Figure [7| 



D. Higher Orders 

In the expansion of the Pauli spin matrices of the last section we have neglected all generated operators with more 
than one bosonic operator. In order to confirm that the expansion of the Pauli spin matrices in normal ordered 
bosonic operators is indeed systematic we will now upgrade our expansion and also include: 

• all generated operators up to two normal ordered bosonic operators 

• all generated operators up to three normal ordered bosonic operators 

In the following, normal ordering shall be defined with respect to the bilinear Hamiltonian of the un-shifted mode, i.e. 
5 = 0. This will simplify matters considerably. Choosing the parameters of the initial Hamiltonian such that (5 = 
for all £, we are still consistent within our normal ordering procedure. 

The first extension, er™ eu '' 2 , includes the following terms, where we introduce the abbreviations x = b + and 
p = b — b^ and where we also confine ourself to the discussion of a z in order to drop one index: 

a new,2 = + ax ^x,+ . ^2 . +i(Ty ^V,+ . xp . 

+ a z 4> z < + : x 2 : +a x ^ x ^ : p 2 : +a z ^ z ^ : p 2 : (62) 
The second extension, er™ eiu ' 3 , consists of the following terms: 

a ne W ,3 = ^1,+ . x 2 . . p 2 . +(Tx{p *,+ ■ x * ■ +ia y (p^+ : x 2 p \ 

+ a z (p z ' + : x 3 : +a x tp x '~ : xp 2 : +ia y <p v ~ : p 3 : +a z ijf~ : xp 2 : (63) 
The resulting flow equations for the upgraded truncation schemes are presented in Appendix IU1 



E. Numerical Results 



We are now set to compare the fixed points of the operator flow obtained from the flow equation approach with 
the exact results. We can also see from the exact solution if the expansion into normal ordered bosonic operators is 
preferable. 

It turns out that the expansion into normal ordered operators is not without obstacles. Especially when the reflection 
symmetry is broken, i.e. eo ^ 0, the final values of the coefficients delicately depend on the initial parameters of the 
Hamiltonian. The reason for this is that the unperturbed states cross when the interaction is switched on and this 
effects the representation of the operator. The effect is enhanced by explicitly breaking certain symmetries. 

Also the comparison of the operator flow with respect to the different versions of the flow equations, discussed in the 
previous section, is troublesome. Since the non-trivial versions for e = are based on different diagonal Hamiltonians 
Hq, a direct comparison of the fixed point parameters is not obvious. 

We therefore limit our investigations to the parameter regime where the reflection symmetry is not broken, i.e. 
eq = 0. If we choose the generator of Version a with ip = and 9 = 0, 5 = for all £, and if we consider the flow of the 
^-component of the Pauli spin matrices, only two parameters h z and x x ' z arc being renormalized. The final values 
h* z = h z (£ = oo) and x x ' z * = X x ' z {^ — °°) are shown for the initial condition Ao/t^o = 0.5 in Figure^ together with 
the results where we also included the flow of bilinear (2. order) and trilincar (3. order) bosonic operators, governed 
by Eqs. - iUl3t and Eqs. itHTH) - (IH27ll . 

The fixed point coefficients h* z and x x ' z * agree with the exact solution unless the initial tunnel-matrix element Ao 
is close to a resonance, i.e. Aq ~ u>o or Aq « 3wcp£. The spike at Aq « 3wq cannot be accounted for by any of the 
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FIG. 8: The fixed point parameters h* = h z (l = do) (left hand side) and x x,z * = X x ' z (^ = 00 ) stemming from the symmetric 
flow equations of Version a for ij) = 0, 9 — 0, Ao/wo = 0.5 and eo = for different orders of truncation of the operator flow as a 
function of Aq. The solid lines resembles the exact result. 



solutions obtained via flow equations. But there is a significant improvement from the second order to the first order 
result close to the resonance at Ao ~ <^o especially in the case of x x ' z * ■ The improvement from third to second order 
in the case of x x,z * is n °t as strong and the one particle parameter h* is almost left unchanged. 

In Figure the results for the fixed point operator a* are shown as they follow from the flow equations of Version 
a with the initial conditions \q/luq — 0.5 and eo = 0. Four parameters g* = g x (£ = 00), /* = f x (£ = 00), 
X v,x * = x v ' x (£ = 00), and x z ' x * = X z ' x (^ = 00 ) are generated during the flow. They show the same deficiencies with 
respect to the exact solution as the results of Figure [H] We want to mention that the constant term /* is indeed 
generated, as can be seen from the exact expansion. 

To investigate the reason for the above discrepancies close to resonances further, we are going to employ the 
numerically exact solution and determine the expansion of the final operator a* = Ua z U^ including up to nine 
bosonic operators. Instead of analyzing the graphs of all 115 coefficients, we will consider the sum of the absolute 
values of the coefficients that belong to the operator class which consists of n bosonic operators (nth-order) . 

The resulting nine graphs are shown in Figure ITU1 As can be seen, the second order still contributes to the fixed 
point operator considerably. Close to resonances even higher orders become important for the operator expansion. 
This explains why the fixed point parameter h* is not sufficiently recovered by the flow equation approach even after 
including all terms up to three bosonic operators into the flow equations. 

In Appendix [DJ the spikes of Figure ^| are related to degeneracies. The formalism thus breaks down at these 
parameter configurations. This is related to the problem that occurs when diagonalizing the Hamiltonian which is, 
strictly speaking, also only possible for non-degenerate states. 

Let us finally comment on sum rules that stem from operator relations which remain invariant under unitary 
transformations. Taking the initial values for the Hamiltonian as in Figure the flow equations of Version a yield 
the exact sum rule (a 2 ) — h 2 + (x x,z ) 2 — 1 at T — for all £ and independent of the initial tunnel-matrix element 
Aq. The sum rule is thus not sensitive to the deviations between the flow equation results and the exact solution, 
which become especially drastic close to resonances, see Fig. |S] We observe the situation that two errors are being 
canceled to yield the desired result. We therefore conclude that the sum rule cannot be a sufficient criterion for the 
quality of the operator flow. On the other hand, one cannot expect that the flow equations yield good results on all 
energy scales. Properties at low energies like the ground-state expectation value of a z shown in Figs. [3] and still 
can be calculated with high precision. The typical deviations at resonances in the operator flow are averaged out. 
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FIG. 9: The parameters g* = g x (t = oo) and f£ = f x (i = oo) (left hand side) as well as x y ' x * — X v ' x (^ = 00 ) an( i 
X z ' x * = X Z,:E (^ = °°) stemming from the symmetric flow equations of Version a for ij) = 0, 9 = 0, Ao/t^o = 0.5 and eo = as 
function of Aq. The solid lines resemble the analytic results. 



IV. CONCLUSIONS 

This work addresses general questions concerning the flow equation approach such as optimization of the final results 
or invariance with respect to the initial Hamiltonian, based on a simple non-trivial model. The model is structurally 
similar to quantum impurity models and since the "bath" only consists of one mode, it is numerically exactly solvable. 
We intended to demonstrate that a systematic improvement of the flow equation approach is possible. In order to 
improve the flow equations one can basically precede according to the following lines: 

1. Most obviously, one can include more interaction terms in the truncation scheme of Hamiltonian and operator. 
This was done for the Hamiltonian flow when employing the form-invariant truncation scheme and a systematic 
improvement was seen. We did not extend the truncation scheme for the canonical generator because it is in 
principle not feasible for more realistic models with an arbitrary number of bosonic modes. For the operator 
flow, the truncation scheme was extended up to third order for a special parameter regime and the results were 
compared with the exact solution on the operator level. Close to resonances, the flow equation results showed 
significant deviations with respect to the exact solution . These deviations were present even in the upgraded 
truncation schemes since high orders of up to nine bosonic operators still carried considerable weight. This is 
connected to the general problem that the flow equation approach breaks down close to degenerate states. 

2. Another way to improve the flow equations is to consider the neglected operators more thoroughly, i.e. to 
introduce a refined decoupling scheme. This was done by introducing a ^-dependent bosonic shift 8 and neglecting 
normal ordered bilinear bosonic operators with respect to this shifted mode. The bosonic shift was deduced 
by formally diagonalizing the truncated Hamiltonian and than decoupling the "system" from the "bath" . The 
decoupling process was not unambiguous and two different approaches were investigated. These were labeled as 
direct and self-consistent evaluation of the system expectation values. The self-consistent approach turned out 
to yield better results on the level of the Hamiltonian flow, the direct approach was preferable on the level of 
the operator flow. 

3. A third possibility to obtain better results is to choose a different basis which the flow is defined on. As an 
example we want to mention the vertex flow introduced by KehreinS*. We investigated the operator flow with 
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FIG. 10: The sum of the absolute values of the coefficients of all operators that consist of n bosonic operators (nth order) 
which compose a*, for Xq/ujo = 0.5 and eo = as function of Aq. 



respect to the distinguished bosonic mode b = b + 5/2. The infinitesimal unitary transformations are then 
equivalent to an active and passive transformation since the coefficients as well as the operator basis b are 
changing during the flow. But the numerical results turned out to be worse than the ones based on the flow 
with respect to the unshifted mode b. We therefore did not include them in the discussion of the present paper. 
We want to emphasis, though, that there remains the possibility to improve the flow equation results along 
these lines. 

It is also pointed out that flow equations are, in general, not invariant with respect to the initial Hamiltonian even 
though the Hamiltonians only differ by a unitary transformation. We concluded that differences are, in general, small 
and if one chooses a form-invariant truncation scheme, the flow equations might not differ at all. But the fact that the 
results depend on the unitary representation of the initial Hamiltonian opens up the possibility to optimize the results 
by introducing an (arbitrary) number of parameters associated with possible unitary transformations and choosing 
them such that certain sum rules are fulfilled best. This strategy has been applied to the spin-boson model with 
external bias, where one parameter - associated with the shift of the bosonic operators - was chosen such that the 
sum rule of a z was optimal for all i£m What had to be left open was how to choose the optimal initial Hamiltonian 
for the evaluation of a specific quantity - a priori. 

The last part of the paper is dedicated to a detailed analysis of the operator flow. Since the flow equations are 
usually designed such that the Hamiltonian is diagonalized best, i.e. that the flow only involves few flow parameters, 
the transformation of the observable is more susceptible to uncontrolled approximations, i.e. higher order interaction 
terms are often neglected merely because they cannot be kept track of. For this reason, the exact operator fixed 
point was evaluated, represented in the basis which was determined by the specific choice of the generator. It turned 
out that the flow equations of the operator should include up to 115 interaction terms in order to adequately coincide 
with the exact operator fixed point on all energy scales. We also pointed out that exact sum rules resulting from the 
flow equations are mostly due to high symmetries of the operator flow, i.e. when only few terms are being generated. 
The assumption that the flow is well approximated if a sum rule holds can thus be misleading as was shown in the 
last section. Nevertheless, the deviations at points of degeneracies of the operator flow with respect to the exact 
solution are unimportant for the low energy properties of the system. This was demonstrated by evaluating the 
ground-state expectation value of o~ z within the most simple, but non-trivial truncation scheme. 
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APPENDIX A: THE INDEPENDENT BOSON MODEL 



We want to give a brief introduction to the flow equation method based on the exactly solvable Independent Boson 
Model. The Hamiltonian of this model is given by 

H = # + ^ = ^ t & + ec t c + Ac t c(& + & t ) . (Al) 

The resemble bosonic, the fermionic operators. They obey the canonical commutation and anti-commutation 
relations respectively. The model can account for some relaxation phenomena and is extensively discussed in the 
textbook by Mahan 2 ^. 

We set e = A 2 /w. Then the Hamiltonian of Eq. IjAip is unitarily equivalent to H = Lob^b + <r z \(b + b*), where a z 
denotes the ^-component of the Pauli spin matrices. This is the Rabi Hamiltonian JJJ with Ao = 0. 



The model is easily solved by the unitary transformation 

£7 = exp(-c t e- (&-&*)) (A2) 

LO 

and we obtain the diagonalized Hamiltonian UHW — ujwb. 

But we want to perform this unitary transformation continuously by introducing a flow parameter £ and a family 
of unitarily equivalent Hamiltonians H(£) = U(£)HU> {£), We also want to look closely at the transformed operator 
c{£) = U(£)cU' (£) and question if an expansion of the operator in a series of unbounded operators, namely (6 — b') n , 
is well-defined. 

The unitary operators U {£) shall be defined by the generator r\ which governs the differential form of a continuous 
unitary transformations as follows: diH = [r),H]. A good choice for the generator has proven to be 77 = [ifo,^], 
which is likely to eliminate the interaction in the limit £ — > oo£. The ^-dependent unitary operator U (£) is 
related to the generator r\ through the differential equation dfU — rjU which can be formally integrated to yield 
U{£) = £exp( J* d£'r] (£')). The operator C denotes the ^-ordering operator, defined in the same way as the more 
familiar time-ordering operator T. In fact, the differential form of the flow equations has got the same structure as 
the Heisenberg equation of motion, but complete formal equivalence is only achieved for explicitly time-dependent 
Hamiltonians, since the generator rj is explicitly ^-dependent. 



For the independent boson model the canonical generator reads r\ = — uj\c^c{b — b*) and we readily obtain 

[r),H] =-cj 2 Ac t c(fe + 6 t )-2cjA 2 c t c . (A3) 

The following flow equations 

di\ = -uj 2 \ , d e e = -2uj\ 2 (A4) 

are integrated to yield X(£) = Xexp(—uj 2 £) and e(£) — — exp(— 2uj 2 £). Since [r)(£), ?](£')] = 0, the ^-ordering operator 
L becomes trivial and we obtain for the ^-dependent unitary operator 

U(£) =exp(-c t c-(l-e^ 2f )(6-fe t )) . (A5) 

LO 

From Eq. I|A5|I we can obtain the unique unitary operator for £ — *■ 00 which diagonalizes H and which was already 
given in Eq. (|A2|) . 

Given U(£) one can determine the flow of the operator c(£) directly: 

c{£) = U(£)cUU£) = cexp(^^(& - &+)) (A6) 

LO 

= cexp(-- exp ^6t) exp (^2 & ) (A7) 

2 \ LO J LO LO 

to J n\ 



v ' n=Q v 
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where we introduced SX(£) = A(l — e~ u e ) and denned normal ordering, denoted by : ... :, by writing the creation 
operator left from the annihilation operator. This definition of normal ordering resembles a special case of the general 
definition given in Appendix A and is valid at T = 0. But from now on the general definition will be used. 

We now apply the continuous transformation to the operator c using the differential form die = [77, c]. The flow 
equations generate the infinite series c(£) — cJ2 n =o ln{£)(p — with c^7„+i = ujX(£)^ n . Together with the initial 
condition 70 = 1, j n = for n > 1, this set of differential equations can be solved to yield j n — -^r( )"■ The flow 
equation result thus coincides with the non-normal ordered form of c(£) in Eq. I|A6() if one expands the exponential 
function into a Taylor-series. 

At first sight there is no distinguished expansion of c{£) in bosonic operators since its generation depends on 
rj. In order to discuss a different scheme, we now define c{£) by a series of normal ordered operators, i.e. c{£) = 
c En=o^W : — ^ )™ ■■ W e obtain the following flow equations 

dan+i = u\(t)(7n - (n + 2)7^+2) , (A9) 

where we used the formula (see Appendix IE)l 

(6 - 6 f ) : (6 - tf) n :=: {b - 6 t )" +1 : +n((b - tff) : (b - 6 t )"- 1 : (Af 0) 

at T = 0, i.e. ((b — b^) 2 ) = —1 with (...) denoting the bosonic ground-state expectation value. Taking the same initial 
conditions as in the case of the non-normal ordered expansion, we see that the normal ordered expansion in Eq. (|A8|) 
solves the set of differential equations (|A9|) . i.e. 7„ = exp(— ^(SX(£)/uj) 2 )^(^p-) n . 

This is a remarkable result. Whereas the non-normal ordered expansion of c{£) reproduces the perturbative result 
in the coupling 8X for each coefficient 7„, the normal ordered expansion yields coefficients 7„, which contain all 
powers of 5X. Especially in view of later approximations, the normal ordered version will then be more preferable, 
since it is likely to go beyond a perturbative description. 

After having recovered the correct flow of the observable via the flow equation approach, we would like to investigate 
the "stability" of the infinite expansion of c(£) in unbounded operators. For this purpose, we consider the Green 
function G(t) = —i(Tc(t)c') and the spectral function A(uj) = — \vaG{ui)J'K with the time ordering operator T, the 
Fourier transform G(ui) = J dte lult G{t) and (...) denoting the ground-state expectation value with respect to H . 
With A e A/w we obtain^ 

G{t) = ^e(£)cxp(-A 2 (l-e-^ t )) , (Afl) 
4(ii) = e 4 TA 2 "^(ii-wj) . (A12) 



n< 

n=0 

The spectral function A(Q) thus exhibits the polaronic shift e p = —\ 2 /u> for n — and an equidistant satellite 
structure separated by the oscillator frequency u> with exponentially decreasing weight. 
Using flow equations, the Green function is best expressed as 

G{t) = -iQ{t){e iH ^°°>c\l = oo)e- lH< - i=co ^c(£ = 00)) , (A13) 

because then the time evolution of the fermionic and bosonic operator is that of free ones. 

In order to recover the exact result, we first use the normal ordered expansion of c(£). With D{t) = b(t) — b^(t), 
where the time evolution is given by the Heisenberg representation with H(£ = 00) = uitfb, the Green function reads: 

G{t) = ~ t e(t)e-" (c(t) £ - : D n (t) : C t ^ —(-1)™ : D m (0) :) (A14) 

n—0 m— 

= -iQ(t)e- x Y, - — (^r(--D n (t),.D m (Q):) (A15) 
L — ' n! ml 

n,m— 

~2 A™ A™ 1 

= -iS{t)e- x - V n \5 ntm e- mM (A16) 

n,m— 

To get from Eq. (|A15ft to Eq. (|A16|) we used the following formula ( Appendix lE|l : 

^fi-fi+r-Cfi-fttr^rexp^-ftt) 2 )^^)^!^^.^: , (A17) 
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with ((b — b^) 2 ) = —1 and (: (b — b') n :) = at T = 0. Summing up the series in Eq. i|A16(l indeed yields the exact 
result given in Eq. I jAllj l. 

In order to show that also the non-normal ordered expansion of c(t) leads to the correct result, we have to normal 
order this expansion. For this we need the following formula (Appendix lE|l : 

(b-^T= J2 fci 2 ^n-2fcV ((b ' bt)2)fc:(b ' bt) "' 2fc: (A18) 
fc=0 ' v '' 

Considering for the moment only the first (N + 1) even powers of (b — b'), we obtain 

N \2n An 7 2 k n k \2(n-k) 

n=0 n=0fc=0 v ' 

-Esff^-"'*^^ ■ (A20 » 

m=0 fc=0 

where we introduced G = ((b — b') 2 ), (...) denoting the canonical ensemble average over a free bosonic system . The 
summation of the first (N + 1) odd powers of (b — b') yields 

N \2n+l N \2m+l N ~ m \2k rik 



n=0 v ; m=Q y > k=0 

In the limit N — > oo we obtain 



n\ e — ; n 

n=0 n=0 



which is an extension of the previous normal ordering of Eq. I|A8(1 to finite temperatures, since G — — (1 + n), 
n = (e^ u — 1) _1 being the Bose factor. This shows that both expansions of c(£) are equivalent. 

To complete the discussion we will now verify that the anti-commutation relation {c(£),c^(£)} = 1 holds for all £. 
To show this we will employ the non-normal ordered expansion. This yields 

Wlc\l)} = \ 7«7»'((-l)" + (-!)*')(& -& f r + "' (A23) 

n,n'— 
2n 

^^(-ljVa^^^i ■ (A24) 

n=0 fc=0 

Summarizing, the series expansion of an operator into bosonic operators yields consistent results. This is no trivial 
result since expanding the bounded operator c into unbounded operators (b — 6')" might have led to inconsistencies. 
Further, it has to be born in mind that the initial operator of the operator flow is resembled by c(£ — 0) = c <g) 1b, 
with 1b being the unity operator of the bosonic Hilbert space. One consequence then is that the trace of the initial 
operator is unbounded and thus not defined. 

As a second result, we want to mention that both expansions, normal ordered and non-normal ordered, are equivalent 
if no approximations are involved. Nevertheless, the operator expansion into normal ordered operators seems to be a 
distinguished expansion since it resembles a non-perturbative approach including the Debye- Waller factor. 

APPENDIX B: THE CONSTANT TERM IN THE EXPANSION OF a x AND a z 

In this appendix we want to comment on the constant term /j appearing in the expansion of the Pauli spin matrices 
o~ x and This term seems to contradict the theorem of the invariance of the trace under unitary transformations. 
But since the trace of Ui{l = 0) = <7j <g> 1, 1 being the identity of the bosonic Hilbert space, does not exist and since we 
also expand the Pauli spin matrices in a series of unbounded operators the above mentioned theorem does not hold 
anymore. To make sure that the constant term is indeed physical, one can truncate the Hilbert space by introducing 
the "bosonic" operator 



>N 



= 6 v /(l-6tfe/JV) , (Bl) 
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with N being a positive integer. The truncated Hilbert space is now only spanned by N vectors \v) = (ft) v /\v\\G) 
with v = 0...N — 1 and 6|0) = 0. For N — > oo we recover the bosonic Hilbert space. The above theorem is guaranteed 
due to the new, non-canonical commutation relation [6/v, frjy] = 1 — (1 + 2b' f b)/N which obeys the cyclic invariance of 
the trace: 

N-l „ 

tr([6 W ,6] v ]) = E( 1 -^^) = ( B2 ) 

The flow equations now have to be extended to include the flow of the operator Wb that appears in the commutator 
relation and that scales as 1 /N. The constant term fi appears nevertheless and is governed by the same differential 
equation as N — > oo. Both terms together, the constant term fi and the bosonic bilinear b^b, make sure that no trace 
is generated during the flow. 



APPENDIX C: UPGRADED FLOW EQUATIONS FOR a z 

In this Appendix, we will set up the flow equations for the Pauli matrix a z including higher orders in the bosonic 
operators. Since the basic objects of our expansion are normal ordered operators we will first give some (anti- 
commutation relations which are helpful to evaluate the commutator [77, a z ] (see also Appendix |E|) : 



[ar, : p n x m 
{x, : p n x m : 



:] = — 2n : p 
} = 2:p 



n—l^m 



n m+1 



{p, : p n x m :} = 2:p 



n+1 x m 



[p, : p n x m :] = 2m: p n x m ~ 1 
In 

2n : p n - l x m : 1„ 



-2m : p n x m - 1 



(CI) 
(C2) 
(C3) 



The commutator of two tensor products of the fermionic and bosonic Hilbert space can be written as [oB, o'B'] = 
{o, o'}[B, B']/2+ [o, o']{B, B'}/2 which is a useful identity if the anti-commutator {o, o'} vanishes (o, o 1 £ H. e , B, B 1 G 
H b ). 

We will not present the commutator [r/, a z \ explicitly but only account for the additional terms that appear with 
respect to the previous flow equations. For the first extension cr™ 6 '"' 2 , this yields: 



dig z = 


... + 2r j x X 1 , d e h z = ...+ 


d eX x = 


...-2rf^ y \ n - V^'+l„ 


d eX y = 




di X z = 


... + 2if + V^'+ln 


di X X = 


4?f - W> y + 4rft/j x ' + 


d^ x ' + = 


-2r? X z ~ 2if' v i> z ' + 


deipy = 


2r, z X x - 2r, x X z 


d t r- + = 


2rj y X x + 2r)°> y i] x ' + 




2t] z X v - 2r] ' y ip z <- 


= 


-2 V x x y + 2if< v ip x >- 


to the previous flow equations coming 




, d iX y = - - 4»jV" . dex 




= ... - 2r] z <p y > + l n - 6?7V' + 


d e ip v 


= ... - ir] z lp x 'l n + 4r] x (p z ~ 




= ... + 2ry V + ln + 6??V + 


d e tp x '~ 


= ... -6r) z <p y -l n -4r] y ip z '~ 




= ... + 6r 1 x <p y <-l n + 4r ] y <p x '- 



T new,3 



•ead: 



- 4f? z i/> 1 ' 



(C4) 
(C5) 
(C6) 
(C7) 
(C8) 
(C9) 
(CIO) 
(CH) 
(C12) 
(C13) 



(C14) 
(C15) 
(C16) 
(C17) 
(C18) 
(C19) 
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The flow equations for the new parameters of a™ ew ' yield: 







-f- zvj Lp -f- 077 9? 


[\jZU ) 


oetp 


— *V <P 




(GzTJ 




— —Arj%p 


— z?y '"9? ' 


fr ir > r >\ 

[VjZZ ) 


n y -4- 


= 2rj tp ' 


— 277 V 


(C23) 


d e ^ + 


= 2t7^V z ' + 


+ 2?7 ' V + 


(C24) 


d e ip x ~ 


= 2?f V" - 


2r^>- - 2?7 'V~ 


(C25) 




= 2rf V*' 4 " 


- 2t ? ;i; V z ' + 


(C26) 




= -2rj x iP v 


+ 277^ - + 2r] ' y ip x '- 


(C27) 



APPENDIX D: RABI MODEL IN PERTURBATION THEORY 

In this appendix we will treat the Rabi Hamiltonian in perturbation theory. We want to start from the exactly 
solvable Jaynes-Cummings Hamiltonian which is obtained from the symmetric Rabi Hamiltonian with no bias by 
applying the rotating wave approximation 19 . This approximation neglects coupling or transition terms which are 
energetically unlikely. 

It is useful to write the Hamiltonian in a basis where a x is diagonal. The Rabi Hamiltonian shall thus be given by 
H = ^2 e i c \°i + uo^b + \bc\co + A& f 4ci + X'tfc\c + A'&cjci . (Dl) 

i=0,l 

The operators cj^ and 6^) obey the canonical anti-commutation and commutation relations respectively. We identify 
the Rabi Hamiltonian given in Eq. by setting e\ — eq = Aq and A = A' = 2Ao and the zero external bias. 

The Jaynes-Cummings Hamiltonian is obtained by setting A' = in Eq. (|D 1|) . We want to treat the so called 
off-shell processes, characterized by A', within a systematic perturbation approach. One way to do so is to consider 
the Hamiltonian in the basis {|0; 2n) |1; 2n + 1)} and {|0;2n+ 1), |l;2n)} where the first quantum number resembles 
the fermionic state and the second quantum number the bosonic state. Since the Hamiltonian is symmetric with 
respect to parity the two sets decouple and in the following we will only consider the first set. 

In the above basis, the Hamiltonian is tridiagonal and we define the n-dependent matrices 



noar„\ - (ti + 2 «M) v2n + 1A 
U [n) ~ {V2n~+TX D° s (n + 1] 



D° n {n + l) 



e + (2n+l)cj \j2n + 2A' 



(D2) 



y/2n + 2A' D on (n + 

The determinants can formally be evaluated to yield 

detD on (n) = (ei + 2nuj )detD oS (n + 1) - (2n + l)A 2 det£> on (n + 1) (D3) 
detD oS (n + 1) = (e + (2n + l)w )detL> on (n + 1) - (2n + 2) X ' 2 det D oS (n + 2) . 

The matrix D on (0) resembles the representation of the Rabi Hamiltonian in the above basis. To determine the 
eigenvalue a of the matrix up to 0(X' ) we iterate Eq. (|D3Jl starting with D° n (0): 

det(£ on (0) - fj,) -> [( £l - M )(e + u> - fi) - A 2 ] 

x [(ei + 2lu - u){e + 3uj - p) - 3A 2 ] det(D on (2) - //) (D4) 
-(ei - u)2X' 2 (e + 3w - Ai)det(£> on (2) - /x) = 

For the eigenvalues we make the ansatz /i = //(°) + A /2 /i (1) - There is no linear term in A' since the spectrum of H 
may not depend on the phase of the coupling constant. 

We now order the eigenvalues as follows: The lowest eigenvalues of order 0(A'°), Mo°±i are determined by setting 
the first factor on the right hand side of Eq. (|D4(I zero. We obtain the well-known Jaynes-Cummings result /ig°]_ = 
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€q + ujq — (A =F i?o)/2 with the detuning A = (ex — eo) — w o and the zeroth Rabi frequency Rq = A 2 + 4A 2 . The first 
correction to ^q°j_ then yields 

(i) _ 1 Acjq ± R u; - A 2 

^ ± ^Ti?o2 Wo 2 T i?o^-A 2 ' (U5) 

The result agrees with the perturbative result in the limit A = A' <C A. 

Generally, setting the nth factor of the first line on the right hand side of Eq. I|D4|I zero the nth eigenvalues yield 
a4°± = e o + (2n + l)uj a - (A =F Rn)/2 with i? 2 = A 2 + 4(2n + 1)A 2 . The first correction to /4°± is given by 



Aw ± i?«^o - A 2 Aw T -Rn^o - A 2 

(" + 1 ) o. .2 - 5 ~ 7Z I TTTa + "7 



2w 2 =F fi«wo - (n + 1)A 2 2w 2 ± i?„w - (n - 1)A 2 



(D6) 



The perturbative approach breaks down when degenerated states are involved. This is indicated by the poles 
in the energy corrections (J^lt- Setting the denominator of zero, one obtains for the tunnel-matrix element 

Ao = ujq + \J (2cJq — A 2 ) 2 — 4wqA 2 /wo- Inserting the parameters used for Fig. [SI we obtain Ao « 2.87. This value 
approximately agrees with the value of Ao where the second spike of h z in Figure [S] is seen. 

APPENDIX E: NORMAL ORDERING 

In this appendix we want to summarize basic relations concerning normal ordering. This summary is based on 
unpublished notes by Wegner of the year 2000 in which he presents a general formalism for normal ordering of 
classical and quantum fields with respect to a bilinear Hamiltonian 26 . We will restrain ourselves to the normal 
ordering of bosonic quantum fields. 

Let bk be any linear combination of Bose creation and annihilation operators. The matrix G shall describe the 
correlations of the operators b for a Hamiltonian H bilinear in the creation and annihilation operators: (bkbi) = Gm- 
The commutator is given by [bk, h] = Gki — Gik- Normal ordering of an operator A with respect to the Hamiltonian 
H is now defined by: 

: 1 : = 1 (El) 
: aA{b) + l3B{b) : = a : A(b) : +/? : B(b) : (E2) 

b k : A(b) : =: b k A(b) : + £ G M : : (E3) 

The product of m operators bk t with i = l..m is now obtained by iterating the third equation. One obtains 

b kl h 2 ...b km =: (b kl +J2Gkuh^Kh 2 +J2 G ^h^)---h m : , (E4) 

h h 

which can also be written as 

b kl bk 2 -b km =■ exp(^G M , e/ ght )b kl bk 2 -b km ■ ■ (E5) 
ki db k db i 

This is Wick's first theorem^. The superscripts le ^ f and right indicate that we always pick a pair of factors b and 
perform the derivative d/dbk on the left factor and the derivative d/dbi on the right factor, so that the factor Gki 
depends on the sequence of the operators. 
Similarly one obtains 

■ b kl b k2 -bk m ■■= exp(- °ki dh left db ri B ht )bhib k2 ...b km . (E6) 
The formula for the product of two normal ordered operators is given by 

: A(b) :: B{b) :=: exp(]T G kl ^j— ) A(b)B(a) : \ a=b . (E7) 

kl k 0,1 
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This is Wicks's second theorem. 

One can now show that under normal ordering the commutative law holds: : ABCD :=: ACBD :. This is rule C 
of Wick. 
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